SMML estimators for 1-dimensional continuous data 



James G. Dowty 

"(N" 

3 '• December 21, 2012 

(N 
o 

Q ' Abstract 

o 

A method is given for calculating the strict minimum message length (SMML) estimator 
for 1-dimensional exponential families with continuous sufficient statistics. A set of n 
equations are found that the n cut-points of the SMML estimator must satisfy. These 
equations can be solved using Newton's method and this approach is used to produce new 
^ ! results and to replicate results that C. S. Wallace obtained using his boundary rules for 

the SMML estimator. A rigorous proof is also given that, despite being composed of step 
functions, the posterior probability corresponding to the SMML estimator is a continuous 
function of the data. 

^ . 

^ ■ 1 Introduction 



The minimum message length (MML) principle [3] is an information theoretic criterion that 
links data compression with statistical inference [3] . It has a number of useful properties 
■ and it has close connections with Kolmogorov complexity [S]. Using the MML principle 

to construct estimators is known to be NP-hard in general [T] so it is common to use 
approximations in practice The term 'strict minimum message length' (SMML) is 
- r-j _ used to distinguish the exact MML criterion from these approximations. 

rS ' The only known algorithm for calculating an SMML estimator is Farr's algorithm [1] 

I which applies to data taking values in a finite set which is (in some sense) 1-dimensional. 

" ■ ■ For 1-dimensional continuous data, certain rules of thumb called boundary rules can some- 

times be used for calculating the SMML estimator [3]. However, these rules were derived 
from a heuristic criterion and are not in general satisfied by the SMML estimator. There- 
fore the calculation of the SMML estimator, even in the simple case of 1-dimensional 
continuous data, is an open problem. 

This paper gives a method for calculating the SMML estimator for a 1-dimensional 
exponential family of statistical models with a continuous sufficient statistic. Section [5] 
recalls the relevant definitions and fixes our notation. Our main results appear in Section 
[3J where we give equations that the cut-points of the SMML estimator must satisfy, show 
how to solve these equations with Newton's method and prove a previously unknown fact 
about the SMML estimator. These results are based on certain technical lemmas whose 
proofs are deferred to Appendix [Al We then apply the results of Section |3] to examples (in 
Sections m and El) before addressing some numerical issues (in Section |6|). Section [7] states 
our main conclusions and discusses some ideas for further research. 
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2 The SMML estimator 



In order to define our notation, tliis section briefly recalls the definition of the SMML 
estimator for a 1-dimensional exponential family of statistical models with a continuous 
sufficient statistic. 

Let the exponential family have support X and natural parameter space Q and assume 
that both are open, connected subsets of E. For each 9 G Q, let f{x\9) be the probability 
density function (PDF) on X given by 

f{x\e) = exp{xe - i^{9))h{x) (1) 

for any x € X, where and h are given functions with h strictly positive everywhere on 
X. If TT is a Bayesian prior on O then we define the marginal PDF r to be given by 

r{x) f ^{d)f{x\9)d9 

for any x d X^ and r{x) — elsewhere. We make the technical assumption that the first 
moment of r exists. 

For the 1-dimensional case considered above, the SMML estimator with n cut-points is 
defined as follows [3]. Suppose we are given an integer n > 1 and real numbers ai < . . . < 
an in X (the cut-points) as well as 9q^ . . . , 9ji € (the assertions) and qq^qi-, . . . Qn ^ E 
(the coding probabilities for the assertions) so that 1 = qq + qi + . . . + qn and each qi > 0. 
Then for each i = 0,1, ... ,n, define Ui to be the interval Ui = [0^,0^+1) where ao and 
a„+i are the boundaries of X, e.g. if A" — M then ag = — cxd and a„+i — 00. Let 9 and q 
be the step functions given by 9{x) =' 9i and q{x) qi where i is the unique integer for 
which X G Ui. If we discretize the data space X to a lattice then there is a 2-part coding 
of the data which has expected length 

h = -E[\ogiqiX)fiX\9iXm (2) 

plus a constant which only depends on the width of the lattice, where X is a random 
variable with PDF r, written X ^ r. Then an SMML estimator with n cut-points is a 
function 9{x) which minimizes Ii out of all estimators of this form. 

This minimality condition can be used to solve for the assertions and the coding prob- 
abilities in terms of the cut-points. Let : 8 — >■ M be the function /i(0) =^ J^xf{x\9)dx 
which relates the natural parametrization of the exponential family to the expectation 
parametrization. By a standard result for exponential families (e.g. see Theorem 2.2.1 
of P]), tp is infinitely differentiable, fi = ip' and fi has an infinitely differentiable inverse. 
Then it is not too hard to show (see R2 and R3 on pages 155-156 and 168-169 of [3]), for 
each i = 0,1, . . . ,n, that 

qi — / r{x)dx (3) 

JUi 

and 

9^ = (^^ J xr{x)dx^ . (4) 

So ^ says that qi is the mass of Ui and ^ says that the centre of mass of Ui is the 
expectation parameter corresponding to 9i . 

Note that an SMML estimator with n cut-points might not exist or might not be 
unique in general. However, we will often refer to 'the' SMML estimator when discussing 
this estimator informally. 
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3 Constructing the SMML estimator 



This section describes our construction of the SMML estimator. This construction is given 
in terms of the natural parametrization of the exponential family but this determines the 
SMML estimator in general since this estimator transforms simply under reparametriza- 
tion. 

Using ([3]), (HJ and the fact that Ui — [0^,0^+1), we can consider qi and 9i to be functions 

of the cut-points a (ai, . . . , a„) G M". Then /i becomes a function solely of a and each 
SMML estimator with n cut-points corresponds to a value of a £ R" which minimizes 
this function Ii{a). But r is continuous so qi{a) and 0i{a) are continuously differentiable 
functions of a, hence so is /i(a) by ^ below. Then since /i(a) is defined on the open 
subset 

R^ = {a e A-" I ai < . . . < a„} 

of R", its gradient vanishes at its minimum (if a minimum exists, i.e. if an SMML estimator 
with n cut-points exists). For each j = 1, . . . , n we therefore have an equation 

|^i = (5) 

which is satisfied at any a € R" corresponding to an SMML estimator. These n equations 
can then be used to solve for the n unknowns oi, . . . , a„, giving the corresponding SMML 
estimator by @ and (H)). 

The next lemma therefore calculates the partial derivatives which appear in ([5|). 

Lemma 1. Let qi and 9i be the functions of a ^ R" given by (0) and ^ and let C be the 
constant — r{x) log h(x)dx. Then 

n 

h{a)^C~Y. * Oog'?^ + kl^m - m)) (6) 

and, for j = l,...,n, 

Proof. See Appendix [XI □ 

Note that the numerator and denominator in the logarithm of ([7]) are, respectively, the 
limits of q{x) f {x\d (x)) as x approaches aj from above and below. Therefore ([5]) is exactly 
the condition which ensures that q{x)f{x\6{x)) is a continuous function of x at a; = Oj. So 
even though q{x) and 9{x) are step functions, we have proved the following. 

Corollary 2. For the SMML estimator, q{x)f{x\9{x)) is a continuous function of x. 

Now, let G : R" — > R" be the function whose j*'* co-ordinate is given by 




G,-(a)=log| -^Vr . I (8) 
for any a G R" . By Lemma [1] 

^^=r{a,)GM, 

so since r{aj) is never zero, solving the system of equations ([S]) is equivalent to the simpler 
and numerically better-behaved problem of finding the zeroes of the function G : R" R". 
We will use Newton's method to find the zeroes of G so the next lemma calculates the 
Jacobian matrix of G and shows that it is sparse. 
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Lemma 3. For j,k = 1, . . . ,n, 

tf\j~k\ > 1. 



dak 



dG 



9G 



9G 



Proof. See Appendix El □ 

Note that n'{9) is the variance of the distribution ^ for any 6 & Q (e.g. see Theorem 
2.2.1 of [2]). 

Remark 1. A global minimum of Ii : K" — )■ R zs ('</ie sei o/ cut-points of) an SMML 
estimator, but solutions to the system of equations (O are only critical points of h . We 
can use Lemma\M to check if a solution to ^ is a local minimum, but these might not be 
global minima. 

Remark 2. // an SMML estimator with n cut-points a*-"-* S M" exists for each n then 
/i(a''"-') is a non-increasing function ofn. To see this, note that /i(a'^"^) < /i(a) for every 
a e R" since a^") is a global minimum of Ii : R" — > R. Also, there exist a e R" with Ii{a) 
arbitrarily close to /i(a("^^)), e.g. take a to be a^"~^' but with an extra cut-point close to 
one of the cut-points of a^"^^^^ and use Therefore /i(a'^"^) < /i(a) = /i(a^"~^)) + e 

for every e > 0, so /i < /i(a("-i)). 

Remark 3. There is a one-to-one map between the set of possible cut-points R" and the 
set of all p £ R" with < pi < . . . < pn < 1, given by pi = R{ai) where R{x) = r{^)d£_ 
is the marginal cumulative distribution function. So we can consider a and hence Ii to 
be a function of p alone, in which case — Gj by ^ and the chain rule, .since the 
Jacobian of the transformation a t-^ p is the diagonal matrix with entries r{ai), . . . , r(a„). 
Parameterizing the cut-points in terms of p has several advantages (e.g. qi is given by the 
simple formula qi — Pi+i — Pi ), but we will not pursue this parametrization here. 



4 Normal data with known variance and a normal 
prior 

We now apply the work of the previous section to a simple case. Each set of cut-points a 
in this section gives a local minimum of Ii{a) but not necessarily a global minimum, so 
we will refer to these as 'likely SMML estimators' to indicate that they are likely but not 
guaranteed to be SMML estimators (see Remark [1]). 

Let A' = = R and choose a normal prior on O with variance , i.e. 9 ^ N{Q,a^). 
Let the data X given 9 be normally distributed with mean 9 and variance 1, i.e. (^|^) ~ 
N{9,1). For example, if Yi,...,F„j are independent and all distributed according to 
N{9^^,a^), where a is known, then (^16*) == is a minimal sufficient statistic for 

Fi,.. and (X|0)~7V(0,1). 

The PDF of X given 9 is of the form (P) with ^P{9) = ^9^ and h{x) = -i=e-"='/2. As 
noted earlier, ^{9) = i}'{9), so ji is the identity map. Also, it is not hard to show that the 
data X (not conditioned on 9) is distributed as X ^ N{0, f3'^) where (3 — y/l + a^, so 
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0.1753120750 


1.9203 


5.9797 
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8840 


17.5442 


27.1130 
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0.1753120750 


1.9203 


5.9797 


10 


8840 


17.5442 


27.1130 


41.1964 
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0.1753120750 


1.9203 


5.9797 


10 


8840 


17.5442 


27.1130 


41.1964 
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0.1753120750 


1.9203 


5.9797 
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8840 


17.5442 


27.1130 


41.1964 


62.1447 
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0.1753120750 


1.9203 


5.9797 


10 


8840 


17.5442 


27.1130 


41.1964 


62.1447 
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0.1753120750 


1.9203 


5.9797 


10 


8840 


17.5442 


27.1130 


41.1964 


62.1447 


93.4500 


16 


0.1753120750 


1.9203 


5.9797 


10 


8840 


17.5442 


27.1130 


41.1964 


62.1447 


93.4500 



Table 1: For various numbers n of cut-points, the difference Ii — Iq in expected code- lengths 
of the one- and two-part codes as weh as the non- negative cut-points bi, . . . ,bk of the hkely 
SMML estimator. 



Table [T] gives, for various numbers n of cut-points, the non- negative cut-points bi, . . . ,bk 
of the likely SMML estimator when a — 2. The bottom line of Table[T](n = 16) corresponds 
to the 'exact SMML' column in Table 3.2 on page 176 of Wallace [Sj and it agrees with 
this column except for Wallace's last entry, which he says is 'not correct'. 

Wallace generated his results using an unspecified iterative procedure which combined 
his boundary rules and ([4]), even though he says these are 'incompatible'. Due to this 
incompatibility, it is maybe not surprising that the boundary rules are not satished for the 
likely SMML estimators given in Table [1] though this makes the close agreement between 
his results and ours even more surprising. It is not clear what connections exist between 
the system of equations (O and Wallace's boundary rules, but there does not seem to be 
a simple connection. 

The SMML estimator seems to be unique and symmetric about when n is 1 or 
3 or n is even, so Table [U determines the likely SMML estimator in these cases, e.g. 
a — {~bk, ■ ■ ■ , —bi,bi, . . . ,bk) if n = 2k. For odd n > 5 there are two likely SMML 
estimators, e.g. if n = 5 the two estimators have cut-points 

a = (-5.9978, -1.9362, 1.9044, 5.9619, 10.8610) 

and 

a = (-10.8610, -5.9619,-1.9044, 1.9362, 5.9978). 

For odd n > 7, each negative cut-point is minus one of the positive cut-points (to four 
decimal places), e.g. when n = 7 the cut-points are 

a = (-10.8840, -5.9797, -1.9203, 1.9203, 5.9797, 10.8840, 17.5442), 

or the negative of this in the reverse order. 

Table [T] also gives the difference /i — Iq in expected code- lengths of the one- and two- 
part codes, where lo = ~ ''(■^) ^^g r{x)dx. Note that increasing the number of cut-points 
beyond n — 6 improves the expected code-length by less than 10~^°, so n = 6 cut-points 
are probably sufhcient for most practical applications. Also note that, to four decimal 
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Figure 1: For the likely SMML estimator with 6 cut-points, the graphs of Dq{x) (dashed line), 
Di[x) (solid line) and r[x) (dotted line) and the corresponding cut-points (vertical lines). 

places, the set of cut-points of each likely SMML estimator with 6 < n < 16 is just a 
subset of the cut-points for the likely SMML estimator when n = 16. This is probably due 
to the fact that more than n = 6 cut-points makes very little difference to Ji and hence 
has little impact on the placement of the existing cut-points. 

Of theoretical interest, there is a local minimum of /i at a = 
(-5.9978,-1.9362,1.9044,5.9619,10.8610,17.5118) which is not a global minimum, 
since /i — /o = 0.1753126143 for this set of cut-points and this is is larger than Ii — Iq for 
the cut-points given in Table [1] for n = 6. So this is a counter-example to idea that all 
local minima of Ji correspond to SMML estimators. 

Figure [1] shows the cut-points and the graphs of Do{x), Di{x) and r{x) corresponding 
to the likely SMML estimator when n = 6, where Do{x) = — r(a:;) log r(x) and Di{x) = 
—r{x)log{q{x)f{x\9{x))) so that Iq — DQ{x)dx and Ii — Di{x)dx. Note that the 
continuity of Di{x), which is guaranteed by Corollary [51 is consistent with this figure. 

5 Exponential data with a gamma prior 

In this section, we apply the results of Section |3] to exponential data with a gamma prior. 

For exponential data, f{x \ 9) is of the form ([T]) if we choose X = [0, oo), = (— oo, 0), 
■il}{9) — — \og{—9) and h{x) = 1 (though note that exponential distributions are usually 
parameterized in terms of the rate —9). Hence the corresponding mean is ^{9) = i/j'{9) = 
— 1/9. Choose a gamma prior for —9 with shape and rate parameters a > and /3 > 
(respectively) so that 

1 a 
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h-h 


ai 


02 


as 


a4 


05 


1 


0.0589128612 


4.49 










2 


0.0579045079 


4.42 


80.55 








3 


0.0579008163 


4.42 


80.17 


1380.63 






4 


0.0579008036 


4.42 


80.17 


1374.66 


23597.96 




5 


0.0579008036 


4.42 


80.17 


1374.64 


23496.46 


403274.23 



Table 2: For various numbers n of cut-points, the difference Ii — Iq in expected code-lengths of 
the one- and two-part codes as well as the cut-points ai, . . . , of the likely SMML estimator. 



Then the data X has marginal PDF 

a I x\ 

i.e. X has a Lomax distribution (equivalently, 1 + X//3 has a Pareto distribution). In 
order to satisfy our technical condition that the first moment of r should exist we need 
to additionally assume that a > 1. Note that for exponential data with an exponential 
prior (a = 1), the expectation defining Ii in 1^ does not in general exist (see (ITUl) ). so the 
SMML estimator is not defined in this case. 

Table [5] gives the cut-points for the likely SMML estimator when a — 2 and (3 — 1. In 
contrast to the normal- normal case of Section 21 for exponential data and a gamma prior 
it seems that the SMML estimator is unique and that all local minima of Ii (a) are global 
minima. 



6 Numerical considerations 

To construct the SMML estimator we might have to consider cut-points ai, . . . , a„ which 
are far outside the likely range of the data, so some of the corresponding values of r{aj) 
and Qj might be extremely small, smaller even than machine precision. This section briefiy 
discusses some simple and effective solutions to the numerical problems that this causes. 



By (H]), the co-ordinate of G : RS. -> M" is given by 



G 



(a) = (logqj + a^Oj - ip{Oj)j - (logq^-i + ajOj^i - -0(^^-1)) 



(9) 



for any a £ R" . For any c € X, let fc{x) = r{x)/r{c), so by ^ we have 



logQj — log ^ J r{x)dx j ~ log ^r(c) J rcix)dx j — log r(c) + log ^ J rc{x)dx j . 

By choosing c appropriately, all terms in the right hand side of this expression can be 
calculated numerically to a high degree of precision (for many functions r{x)). For example, 
with r{x) as in Section|4l fdx) ~ exp(— (a:^ — c'^)/2(3^), so taking c = aj we have 



loggj =-^-^log(2^/32)+log 



1 





"-(x2-a2)" 






dx^ 


2^2 



which is numerically well-behaved even for large aj. Also, by (j4]) we have 
' 1 



\ I frr xr(x)dx\ ( \jj xfcix)dx^ 

xr{x)dx U M^M r . = A^"' ' 



/[/, r{x)dx 



L fc{x)di 
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and the right-hand side is again numerically well-behaved for some choice of c, e.g. with 
r{x) as in Section |4] we could take c = aj if aj > and c = a^+i if Oj < 0. This shows 
that high-precision numerical calculations of Ii and G are possible, even when qj is smaller 
than machine precision. 

We also note that r(oi) and qj only appear in Lemma 12] as ratios of each other. We 
can calculate r{aj+i)/qj by evaluating the right-hand side of 



and this is numerically well-behaved for appropriate c. Other ratios r(a,;)/qj can be cal- 
culated similarly so the Jacobian matrix of G can also be calculated numerically. 

7 Conclusions and extensions 

In the context of 1-dimensional exponential families with continuous sufficient statistics, 
we have found equations that the cut-points of the SMML estimator must satisfy. As 
a corollary, we proved that the posterior probability q{x)f{x\9{x)) corresponding to the 
SMML estimator is a continuous function of x, despite being composed of step functions. 
We also solved these equations for a particular example using Newton's method. Our 
approach is very simple but it solves an outstanding problem in information theory which 
previously could only be attempted with rules of thumb like Wallace's boundary rules. 

Focussing on the case of continuous data allowed us to use calculus to solve the opti- 
mization problem defining the SMML estimator. Restricting to 1-dimensional data allowed 
us to assume a particular form (intervals) for the shape of the regions defining the SMML 
estimator. Therefore our results probably generalize fairly easily to non-exponential fam- 
ilies with 1-dimensional sufficient statistics. It is also possible that they will generalize 
to higher-dimensional continuous data, if the regions which define the SMML estimator 
are assumed to be convex polygons (or any other shapes whose configuration space is a 
manifold) . 

Many questions about SMML estimators for continuous data remain unanswered, even 
in the simple, 1-dimensional case considered here. Does an SMML estimator with a given 
number of cut-points always exist? Is the SMML estimator with a given number of cut- 
points unique for positive data? Does the system of equations ([5]) have a finite number 
of solutions? If the data is restricted to a compact (i.e. finite and closed) interval then is 
there an upper bound to the number of cut-points that an SMML estimator can have? 

An affirmative answer to the last two questions would open the possibility of developing 
a rigorous algorithm to find all SMML estimators with a given number of cut-points (by 
finding all solutions to ([5]) and outputting those with the lowest Ii) and a continuous 
analogue of Farr's algorithm [1] for positive data. 

A Proofs of technical lemmas 

This appendix contains the proofs of our main technical lemmas. We begin with a calcu- 
lation which will be used in both proofs. 

Lemma 4. Let qi and Oi he the functions of the cut-points a given by (0) and Q). Then 
for j = 0, 1, . . . , n and k = 1, . . . , n, 



1j 



lu, r[x)dx 



Su, rc{x)dx 




oq^ 

- — = er(afe) 
oak 



and 
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where 

— 1 if k — i; 

1 ifk = i + l; 

otherwise. 



Proof. Let R be the marginal cumulative distribution function of the data given by R{x) = 
/-co <0d^ for any x € R. Then by ® , = i?(a,+i) - i?(a,;) so = -r(a,), gff^ = 
r(ai-|_i) and = unless k = i,i + 1. 

Now, let = /^^ ^r(C)d^ for any a; G M so that qifi{9i) — Af (a^+i) — M{ai) by 

Differentiating this equation with respect to gives 

09- 

-r{ai)^.{9i) + qi^.'{9i)—-^ = -air{ai) 
da.. 



so by rearranging we have 

39 i r{ai) 



fi{9i) - a, 



da, q,fi'(9,) 

The cases k = i + 1 and k ^ i,i + 1 can be handled similarly. □ 

Proof of Lemma[l[ As in the statement, let C be the constant — r{x) log h{x)dx. From 
([2) we have 



h = - r{x)\og{q{x)f{x\9{x)))dx 
J X 

n „ 

= -E / K^)log(g(x)/(a:|^(a;)))dx 



E / r(a:)log(g,/(x|^,))rfa 



C - V / r(a;) [log <?, + a;^, - V'(^*)]rfa; by (P 

C + Y^^(^q,\ogq, + q,il}{9i) - 9, j xr{x)dx^ (10) 



and © follows by (g]). Then by dTU]) and Lemma H 

5^ ^ X! -g^ y^lt^^SQi + Qii^i^i) - Oi J xr{x)dx 



"(aj) logg^-i - 1 + V'(^j-i) - flj^i-i) + / l^^j-iV^'l^j-i) - J xr{x)dx 

89, 



+r{a,) {\ogq, + 1 - ij{9,) + aj9j^ + ^ [qji^'iO,) - xr{x)dxj . 
Now, xr{x)dx = qjfii9j) by dH) and = by Theorem 2.2.1 of [2 . Therefore the 



bracketed expression multiplying vanishes, as does the expression multiplying 
so we have 

1^ = r{aj) (log qj + aj9j - V'(^j)) - r{aj) (log qj^i + a.j9j^i - ^j{9.j-i)^ III) 



and the lemma follows from ([T]). □ 
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Proof of Lemma\^ By ((T]) and (H)), 

Gj{a) = (loggj + flj^j - - (loggj_i + a^^j-i - V'(^Vi)) (12) 

for any a S M". So = whenever |j — fc| > 1 by Lemma IH 

For the rest of the lemma, just differentiate ([T2l) . use the fact that fi — t// and apply 
Lemma m For example, ii j ^ n then 



dGj ^ I dq, ^ f ^ X de, 1 dq,^, ^ f ^ ^^e,^, 



r 



' (a,^ ti{§,)) li^^ (a^.+i - f^{ej)) by LemmaH 



□ 
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